// ***** This file is automatically generated from LinearTernaryCore.java.jpp

package daikon.inv.ternary.threeScalar;

import daikon.*;
import daikon.inv.*;
import daikon.inv.binary.twoScalar.*;
import daikon.inv.unary.scalar.*;
import plume.*;
import java.io.Serializable;
import java.util.*;
import java.util.logging.Logger;
import java.util.logging.Level;

/**
 * The LinearTernaryCore class is acts as the backend for the
 * invariant (ax + by + cz + d = 0) by processing samples and
 * computing coefficients.  The resulting coefficients a, b, c, and d are
 * mutually relatively prime, and the coefficient a is always p
 **/

// Originally, LinearTernaryCore code was not dealing with degenerate
// linear ternary invariants correct, namely when the plane
// (ax+by+cz+d=0) is parallel to one of the axes, so one of the
// coefficients was 0.  In this case, the invariant can be described
// by a binary invariant (eliminate the 0 coefficient variable).
// LinearTernaryCore's inability to deal with this degenerate form was
// (and still is) masked by the suppression of a linear ternary
// invariant via a binary invariant and a constant (the code threw an
// exception when seeing this case and falsifies the invariant).  As
// more samples are processed, the ternary invariant may "evolve" from
// a degenerate form to a full ternary invariant.  The
// linearintervention method deals with the degenerate form and the
// planarintervention focuses on the full ternary invariant.

@SuppressWarnings("nullness") // ***** TEMPORARY:  def_points is confused (is packed initially with non-nulls at beginning, but is not all non-nulls; and later may not be packed); needs review/refactoring/rewriting
public final class LinearTernaryCore
  implements Serializable, Cloneable
{
  // We are Serializable, so we specify a version to allow changes to
  // method signatures without breaking serialization.  If you add or
  // remove fields, you should change this number to the current date.
  static final long serialVersionUID = 20030822L;

  /** Debug tracer. **/
  static final Logger debug = Logger.getLogger("daikon.inv.ternary.threeScalar.LinearTernaryCore");

  // ax + by + cz + d = 0; first argument is x, second is y, third is z
  // Although the coefficients are guaranteed to be integers, we use doubles to
  // store the values because regression tests have shown that a long is
  // not sufficient
  public double a = 0, b = 0, c = 0, d = 0;
  public double min_a, max_a, min_b, max_b, min_c, max_c, min_d, max_d;
  public double separation = 0;

  // coefficients for the invariant in 3D space
  public double[] coefficients = new double[4];

  // The values of the flags determine indicate whether a line exists
  // with the current samples (NO_LINE) and which variables are
  // constant (constant variables indicate that the 'plane' is
  // parallel to that axis).  For example, XY_CONSTANT means that the
  // so far, the x and y values have all been constant and only z has
  // been varying so so the plane is actually a line.
  // NO_CONSTANT_VARIABLES indicates that the values have not been
  // constant but do for a plane between variables.
  public enum Flag {NO_LINE, XY_CONSTANT, XZ_CONSTANT, YZ_CONSTANT,
                    X_CONSTANT, Y_CONSTANT, Z_CONSTANT, NO_CONSTANT_VARIABLES};

  // The flag indicates whether a line exists in the 3D space instead
  // of a plane If the flag value is not NO_LINE, then a line does
  // currently exist among the values seen and the flag specifies
  // which variables are constant to make the invariant a line rather
  // than a plane.  The value NO_LINE indicates either that not enough
  // values have been seen to make any decision or that a plane
  // exists.
  public Flag line_flag = Flag.NO_LINE;

  // Inner class to store points in for more convenient access.
  public static class Point implements Serializable, Cloneable {
    static final long serialVersionUID = 20050923L;

    public long x;
    public long y;
    public long z;
    public Point () {
      this (0, 0, 0);
    }
    public Point (long x, long y, long z) {
      this.x = x; this.y = y; this.z = z;
    }
    public boolean equals (long x, long y, long z) {
      return ((this.x == x) && (this.y == y) && (this.z == z));
    }
    protected Point clone() throws CloneNotSupportedException {
      return (Point) super.clone();
    }
    public String toString() {
      return ("(" + x + ", " + y + " ," + z + ")");
    }
  }

  // Points that define the plane.  The first three are normally the
  // current definition.  The next is a new point for consideration.
  // However, any slot may be set to null at any time, and points may be
  // duplicated at any time.
  // (The invariant appears to be a bit stronger than that, since there are
  // places that dereference the array unconditionally, assuming that its
  // elements are non-null.  This should be clarified/refactored.)
  public /*@Nullable*/ Point[] def_points = new Point[MINTRIPLES];

  public Invariant wrapper;

  // The number of distinct values (not samples) seen.
  public int values_seen = 0;

  static final int MINTRIPLES = 5;

  public LinearTernaryCore(Invariant wrapper) {
    this.wrapper = wrapper;
  }

  public LinearTernaryCore clone() {
    try {
      LinearTernaryCore result = (LinearTernaryCore) super.clone();
      result.def_points = new Point[MINTRIPLES];
      for (int i=0; i<MINTRIPLES; i++) {
        Point p = def_points[i];
        if (p != null) {
          result.def_points[i] = p.clone();
        }
      }
      return result;
    } catch (CloneNotSupportedException e) {
      throw new Error(); // can't happen
    }
  }

  /**
   * Reorganize our already-seen state as if the variables had shifted
   * order underneath us (rearrangement given by the permutation).
   **/
  public void permute(int[] permutation) {
    assert permutation.length == 3;
    assert ArraysMDE.fn_is_permutation(permutation);
    // Fix a, b, c
    double[] clever = new double[] { a, b, c };
    double[] pclever = new double[3];
    pclever[permutation[0]] = clever[0];
    pclever[permutation[1]] = clever[1];
    pclever[permutation[2]] = clever[2];
    if (isActive()) {
      // Seen enough values, so permuting a, b, c is useful
      if ( a == 0 || b == 0 || c == 0) {
        throw new Error("Active LTs never have 0 coefficients");
        // // Can't handle this form once rotated.  But if we've seen
        // // enough values, yet a, b or c is zero, then this is covered by
        // // a LinearBinary or a constant, so let me just destroy myself.
        // values_seen = Integer.MAX_VALUE;
        // a = b = c = d = 0;
        // if (debug.isLoggable(Level.FINE)) {
        //   debug.fine ("  Ternary invariant destroyed because a, b, or c = 0");
        // }
        // return;
      } else {
        // double d = -1.0 / pclever[2];
        a = pclever[0];
        b = pclever[1];
        c = pclever[2];

        // we still need to guarantee a positive x coefficient
        if (a < 0) {
          a = -a;
          b = -b;
          c = -c;
          d = -d;
        }
      }
    }

    // Fix caches
    {
      long[] temp = new long[3];
      for (int i=0; i<MINTRIPLES; i++) {
        Point p = def_points[i];
        if (p == null)
          continue;
        wrapper.log ("orig def_points["+i+"] = " + p);
        temp[permutation[0]] = p.x;
        temp[permutation[1]] = p.y;
        temp[permutation[2]] = p.z;
        p.x = temp[0];
        p.y = temp[1];
        p.z = temp[2];
        wrapper.log ("permuted def_points["+i+"] = " + p);
      }
    }
    // Assert that caches sync with a,b,c?
    // This would be a good sanity check, but it would be nontrivial
    // because we don't keep track of when a, b, and c are really
    // valid for all the of the points we've cached. If we were
    // falsified and then resurrected, we might have samples that
    // didn't fit even before the permutation. -smcc
  }

  /**
   * Returns whether or not the invariant is currently active.  We become
   * active after MINTRIPLES values have been seen and a line is not calculated.
   * In addition, the coefficients (a, b, c) must all be non-zero, else the
   * invariant would be described by a LinearBinary or constant invariant.
   * Before that, a, b, and c are uninitialized.
   */
  // this is the check used throughout the file as a shortcut because
  // the invariant must have seen enough samples and not be a line for all the
  // other calculations and comparisons done to be sensible, meaning that
  // a full linear ternary invariant exists, not a degenerate one
  public boolean isActive() {
    return (line_flag == Flag.NO_LINE && values_seen >= MINTRIPLES && a != 0 && b != 0 && c != 0);
  }

  /**
   * Sets up the invariant from a LinearBinary invariant and a constant
   * value for the third variable.  Points are taken from the LinearBinary
   * cache and its min_x/y and max_x/y points and combined with the
   * constant value.
   *
   * @return InvariantStatus.NO_CHANGE if the invariant is valid,
   * InvariantStatus.FALSIFIED if one of the points invalidated the
   * LinearTernary invariant
   */
  public InvariantStatus setup (LinearBinary lb, VarInfo con_var,
                                long con_val) {

    if (Debug.logOn())
      wrapper.log ("setup from lb " + lb + " con var " + con_var +
                   " con_val " + con_val);

    int con_index = con_var.varinfo_index;
    int lb_v1_index = lb.ppt.var_infos[0].varinfo_index;
    int lb_v2_index = lb.ppt.var_infos[1].varinfo_index;

    long[] con_vals = new long [lb.core.values_seen+2];
    long[] lb1_vals = new long [lb.core.values_seen+2];
    long[] lb2_vals = new long [lb.core.values_seen+2];

    int mi = lb.core.values_seen;
    for (int i = 0; i < mi; i++) {
      con_vals[i] = con_val;
      lb1_vals[i] = lb.core.x_cache[i];
      lb2_vals[i] = lb.core.y_cache[i];
    }
    if (lb.isActive()) {
      con_vals[mi] = con_val;
      lb1_vals[mi] = lb.core.min_x;
      lb2_vals[mi] = lb.core.min_y;
      con_vals[mi+1] = con_val;
      lb1_vals[mi+1] = lb.core.max_x;
      lb2_vals[mi+1] = lb.core.max_y;
      mi += 2;
    }

    InvariantStatus sts = InvariantStatus.NO_CHANGE;

    for (int i = 0; i < mi; i++ ) {
      if (con_index < lb_v1_index) {
        sts = add_modified (con_vals[i], lb1_vals[i], lb2_vals[i], 1);

      } else if (con_index < lb_v2_index) {
        sts = add_modified (lb1_vals[i], con_vals[i], lb2_vals[i], 1);

      } else {
        sts = add_modified (lb1_vals[i], lb2_vals[i], con_vals[i], 1);

      }
      if (sts != InvariantStatus.NO_CHANGE)
        break;
    }

    if (sts != InvariantStatus.NO_CHANGE) {
      System.out.println ("lb.core.values_seen=" + lb.core.values_seen);
      for (int i = 0; i < mi; i++ )
        System.out.printf ("LTCore: vals %s %s %s%n", con_vals[i], lb1_vals[i],
                           lb2_vals[i]);
      System.out.println ("in inv " + wrapper.format() + " " + wrapper.ppt);
      assert sts == InvariantStatus.NO_CHANGE;
    }

    return (sts);
  }

  /**
   * Sets up the invariant from a OneOf and two constants.
   * Points are taken from the OneOf cache and the constant values.
   *
   * @return InvariantStatus.NO_CHANGE if the invariant is valid, or
   * InvariantStatus.FALSIFIED if one of the points invalidated the
   * LinearTernary invariant
   */

  public InvariantStatus setup (OneOfScalar oo, VarInfo v1, long con1,
                                VarInfo v2, long con2) {

    int oo_index = oo.ppt.var_infos[0].varinfo_index;
    int con1_index = v1.varinfo_index;
    int con2_index = v2.varinfo_index;

    InvariantStatus sts = InvariantStatus.NO_CHANGE;
    for (int i = 0; i < oo.num_elts(); i++ ) {
      if (oo_index < con1_index) {
        sts = add_modified ((((Long)oo.elt(i)).longValue()), con1, con2, 1);
      } else if (oo_index < con2_index) {
        sts = add_modified (con1, (((Long)oo.elt(i)).longValue()), con2, 1);
      } else
        sts = add_modified (con1, con2, (((Long)oo.elt(i)).longValue()), 1);
      if (sts != InvariantStatus.NO_CHANGE)
        break;
    }

    if (Debug.logOn())
      wrapper.log ("setup from OneOf " +  oo + " v1=" + con1 +
                   " v2=" + con2 + " status = " + sts);

    return (sts);
  }

  /**
   * Looks for points that define a plane (ax + by + cz + d = 0).
   * Collects MINTRIPLE points before attempting to define a line
   * through the points.  Once the equation for the line is found,
   * each subsequent point is compared to it.  If the point does not
   * match, recalculates the maximally separated points and attempts
   * to fit the points to the new line.  If the points do not fit the
   * line, attempts to define the plane (to hopefully get at least
   * some spread between the points, so that small errors don't get
   * magnified).  Once the equation for the plane is found, each
   * subsequent point is compared to it.  If the point does not match
   * the point is examined to see if it would is maximally separated
   * when compared to the points originally used to define the plane.
   * If it is, it is used to recalcalulate the coefficients (a, b, c).
   * If those coefficients are relatively close to the original
   * coefficients (within the ratio defined by Global.fuzzy) then the
   * new coefficients are used.
   *
   * @see FuzzyFloat
   */
  public InvariantStatus add_modified(long x, long y, long z, int count) {

    if (Debug.logDetail())
      wrapper.log ("Adding point, x=" + x + " y=" + y + " z=" + z + " to invariant");

    if (values_seen < MINTRIPLES) {
      // We delay computation of a and b until we have seen several triples
      // so that we can compute a and b based on a far-separated triple.  If
      // the points in a triple are nearby, then roundoff errors in the
      // computation of the slope can be non-negligible.

      // skip points we've already seen
      for (int i = 0; i < values_seen; i++)
        if (def_points[i].equals (x, y, z))
          return InvariantStatus.NO_CHANGE;

      def_points[values_seen] = new Point (x, y, z);
      wrapper.log("Add: (" + x + ", " + y + ", " + z + ")");
      values_seen++;
      wrapper.log("Values seen: " + values_seen);

      // check to see if we've seen enough values to create the equation
      if (values_seen == MINTRIPLES) {

        // Try to fit the points to a line first.
        // The points may form a degenerate linear ternary invariants
        // especially if one of the variable is constant
        linearIntervention(def_points);

        // if line can not be formed, try to form a plane
        if (line_flag == Flag.NO_LINE) {
          InvariantStatus stat = planarIntervention(def_points);
          return stat;
        } else {

          // points fit in a line
          return InvariantStatus.NO_CHANGE;
        }
      } else  {

        // still haven't seen enough values
        return InvariantStatus.NO_CHANGE;
      }
    } else {

      //at this point either a line or a plane must have been formed

      // if line already broken, fit to plane equation
      if (line_flag == Flag.NO_LINE) {
        // If the new value doesn't fit the equation
        if (!(Global.fuzzy.eq (-d, a*x+b*y+c*z))) {

          // Try to find small changes that will fit better.
          if (!try_new_equation (x, y, z)) {
            if (Invariant.logOn() || debug.isLoggable(Level.FINE)) {
              wrapper.log (debug, "destroying  (" + wrapper.format()
                           + ") where x=" + x + " y=" + y + " z=" + z
                           + " a=" + a + " b=" + b + " c=" + c
                           + " values_seen=" + values_seen);
            }
            // jiggling the values did not work
            return InvariantStatus.FALSIFIED;
          }
        } else {

          // point fits the current plane equation
          return InvariantStatus.NO_CHANGE;
        }
      } else {

        // try to fit to the current line
        if (!try_points_linear(x, y, z)) {
          //put the crucial point in, so that it can be tested against the new line
          def_points[MINTRIPLES-1] = new Point (x, y, z);

          // try to fit the points onto a line again (not sure if this
          // is really necessary, maybe the fuzzy equals really does
          // make a difference)
          linearIntervention(def_points);
        }

        // new point has broken the line, try to fit to a plane
        if (line_flag == Flag.NO_LINE) {
          InvariantStatus stat = planarIntervention(def_points);
          return stat;
        } else {

          // fuzzy equals does make a difference
          return InvariantStatus.NO_CHANGE;
        }

      }

    }
    // should never get here because values_seen is either < or >= MINTRIPLES
    // return InvariantStatus.NO_CHANGE;
    throw new Error("at end of add_modified");
  }

  /**
   * Attempts to fit the points in def_point to a plane and calculates the
   * coefficients (a, b, c, d) if possible.
   *
   * @param def_point  Array of points.  Must have at least 3 elements.
   *
   * @return the status of the invariant
   *     (whether the plane fitting was successful)
   */
  private InvariantStatus planarIntervention(/*@Nullable*/ Point[] def_point) {

    // make a copy of the array so that the points don't get clobbered
    // because we need to check all the points in def_points against the
    // coefficients calculated

    /*@Nullable*/ Point[] dummy = new /*@Nullable*/ Point[def_point.length];
    for (int i = 0; i < def_point.length; i++)
      dummy[i] = def_point[i];

    // Find the three points with the maximum separation
    maxsep_triples (dummy);
    // Now, first 3 values of dummy are the maximum-separated 3 points.

    // Null points only show up in the cross-checker, not the regression
    // tests, as of 5/25/2010.
    // Sometimes the points are still null.
    if (dummy[0] == null) {
      return InvariantStatus.FALSIFIED;
    }

    // Calculate the coefficients of the equation (a, b, c, d).
    // Requires that the first 3 elements of dummy must be non-null.
    double[] coef = calc_tri_linear (dummy);

    a = coef[0];
    b = coef[1];
    c = coef[2];
    d = coef[3];

    wrapper.log("Calculated tri linear coeff: ("+ a +
                "), b (" + b + "), c (" + c + "), and d (" + d + ")");
    // If one of these coefficients is zero (except for d), this
    // should be a LinearBinary, not a LinearTernary, term.  (It
    // might not show up as LinearBinary because there might not
    // have been enough samples; but a random varying third
    // variable can create enough samples.)  Also, throw out the
    // invariant if any of the coefficients would be not-a-number.
    if ( a == 0|| b == 0 || c == 0 ||
        Double.isNaN(a) || Double.isNaN(b) || Double.isNaN(c)) {

      wrapper.log("problematic coefficient: invariant falsified");
      return InvariantStatus.FALSIFIED;
    }

    // Check all values against a, b, and c.
    if (!wrapper.is_false()) {
      for (int i=0; i<values_seen; i++) {

        // Global.fuzzy.eq works by comparing ratios between numbers,
        // so comparing the sum of everything to 0 won't work.
        //  wrapper.log("point: " + def_point[i].x + " " +
        //  def_point[i].y + " " + def_point[i].z);
        //  wrapper.log("compare: " + -c * def_point[i].z);
        //  wrapper.log( "to: " + (a*def_point[i].x +
        //  b*def_point[i].y + d));

        wrapper.log ("Trying at index " + i + ": "
                     + "0 != "
                     + a + "*" + def_point[i].x
                     + "+" + b + "*" + def_point[i].y + "+" + c
                     + "*" + def_point[i].z
                     + "+" + d);

        if (!(Global.fuzzy.eq (-d, a*def_point[i].x + b*def_point[i].y + c * def_point[i].z))) {

          if (Debug.logOn() || debug.isLoggable(Level.FINE)) {
            wrapper.log (debug, "Destroying at index " + i + ": "
                         + "0 != "
                         + a + "*" + def_point[i].x
                         + "+" + b + "*" + def_point[i].y + "+" + c
                         + "*" + def_point[i].z
                         + "+" + d);
          }

          return InvariantStatus.FALSIFIED;
        }
      }
      if (Debug.logOn())
        wrapper.log ("equation = " + a + "*x " + b + "*y" + c + "*z = "
                     + (-d));

      // Discard the points not used to define the coefficients.
      for (int ii = 3; ii < MINTRIPLES; ii++)
        def_point[ii] = null;
    }

    // if it passes all the checks, then no change to the Invariant
    return InvariantStatus.NO_CHANGE;
  }

  /**
   * Attempts to fit the points in def_point to a line and calculates the
   * the corresponding line classification (line_flag) and coefficients
   * (coefficients[]) if possible.  If not (points in def_point form a plane)
   * resets line_flag to Flag.NO_LINE and coefficients[] elements to -999.
   *
   * @param def_point  Array of points.  Must have at least 2 elements.
   *
   */
  private void linearIntervention(/*@Nullable*/ Point[] def_point) {

    // find the max separation points
    Point[] maxsep_doubles = maxsep_doubles(def_point);

    // maxsep_doubles[0] could be null when enough points have a NaN
    if (maxsep_doubles == null) {
      // line_flag == Flag.NO_LINE so the invariant
      // will try to fit a plane to the points and fail
      return;
    } else {
      // find the equation (coefficients now stored in coefficients[])
      calc_bi_linear(maxsep_doubles[0], maxsep_doubles[1]);

      // This seems to imply that all elements of def_points are non-null.
      // Is that true?
      // check the other points against the coefficients
      for (int i = 0; i < def_point.length; i++) {
        // try to fit the current point to the current line
        boolean ok = try_points_linear(def_point[i].x, def_point[i].y, def_point[i].z);

        if (!ok) {
          line_flag = Flag.NO_LINE;
          break;
        }
      }
    }

  }

  /**
   * Attempts to fit the new point (x,y,z) onto the current line
   * (indicated by line_flag and coefficients[]).  Method should only be called
   * if the current points define a line (line_flag is not Flag.NO_LINE)
   *
   * @param x,y,z  x,y,z components of the point.
   *
   */
  private boolean try_points_linear(double x, double y, double z) {

    switch (line_flag) {
    case XY_CONSTANT:
      return (Global.fuzzy.eq (coefficients[0], x)) && (Global.fuzzy.eq (coefficients[1], y));

    case XZ_CONSTANT:
      return (Global.fuzzy.eq (coefficients[0], x)) && (Global.fuzzy.eq (coefficients[1], z));

    case YZ_CONSTANT:
      return (Global.fuzzy.eq (coefficients[0], y)) && (Global.fuzzy.eq (coefficients[1], z));

    case X_CONSTANT:
      return (Global.fuzzy.eq (coefficients[0], x)) && (Global.fuzzy.eq (z, coefficients[1] * y + coefficients[2]));

    case Y_CONSTANT:
      return (Global.fuzzy.eq (coefficients[0], y)) && (Global.fuzzy.eq (z, coefficients[1] * x + coefficients[2]));

    case Z_CONSTANT:
      return (Global.fuzzy.eq (coefficients[0], z)) && (Global.fuzzy.eq (y, coefficients[1] * x + coefficients[2]));

    case NO_CONSTANT_VARIABLES:
      return (Global.fuzzy.eq (y, coefficients[0] * x + coefficients[1])) && (Global.fuzzy.eq (z, coefficients[2] * x + coefficients[3]));

    // error if the code gets here because that means, we are trying to fit the points to a linear
    // when the points already form a plane
    default:
      throw new RuntimeException("try_points_linear called when points already form a plance (line_flag = NO_LINE)");
    }
  }

  /**
   * Returns the two points that have the maximum separation in pa.
   *
   * @param ps  Array of points.  Must have at least 2 elements.  Can be any
   *            length and can contain nulls (which will be ignored).
   *
   * @return a 2-element array containing the most-separated elements.
   *         Returns null if no pair of elements is NaN-free.
   */
  private Point /*@Nullable*/ [] maxsep_doubles (/*@Nullable*/ Point[] pa) {

    Point p1 = null, p2 = null;

    // max_separation is the separation metric for the most separated
    // pair of points.  We use the sum of the separations as the
    // metric.
    double max_separation = Double.MIN_VALUE;

    for (int i=0; i < pa.length-1; i++) {
      for (int j=i+1; j < pa.length; j++) {
        double separation = separation(pa[i], pa[j]);
        if (separation > max_separation) {
          max_separation = separation;
          p1 = pa[i];
          p2 = pa[j];
        }
      }
    }

    // This can happen if enough values are NaN.
    // if (p1 == null) {
    //   throw new IllegalArgumentException(Arrays.toString(pa));
    // }

    if (Debug.logDetail())
      wrapper.log ("maxsep_doubles = " + p1 + " " + p2);
    if (p1 == null) {
      return null;
    }
    return new Point[] { p1, p2 };
  }

  /**
   * Calculates the coefficients for the line in 3D
   * and sets the coefficients[] by side effect.
   *
   * @param points  2 points to use to calculate the coefficents.
   */
  //
  //  the line can fall under 3 categories
  //
  //  1. parallel to an axis which means that 2 of the following is true
  //  x = A
  //  y = B
  //  z = C
  //
  //  and the two that are true form a line in the 3rd dimension (the one that's
  //  not constant)
  //
  //  2. parallel to the plane of an axis, which means that one of the following
  //  is true
  //  x = A
  //  y = B
  //  z = C
  //
  //  and the line is just a line in the non-constacoefficientble plane projected
  //  onto the third dimension
  //  and the equation for that line is:
  //  x = A [suppose that x is the constant one]
  //  y = Bz + C
  //
  //
  //  3. none of the above and the equation of the line is just:
  //  y = Ax + B
  //  z = Cx + D
  //
  //  Note that A, B, C, and D are coefficients that are stored in coefficients[]

  private void calc_bi_linear(Point p0, Point p1) {

    if ((p0.x == p1.x) && (p0.y == p1.y) && (p0.z != p1.z)) {
      // x and y have constant values
      // line is defined by x = A, y = B
      line_flag = Flag.XY_CONSTANT;
      coefficients[0] = p0.x;
      coefficients[1] = p0.y;
      coefficients[2] = 0;
      coefficients[3] = 0;
      return;
    }
    if ((p0.x == p1.x) && (p0.z == p1.z) && (p0.y != p1.y)) {
      // x and z have constant values
      // line is defined by x = A, z = C
      line_flag = Flag.XZ_CONSTANT;
      coefficients[0] = p0.x;
      coefficients[1] = p0.z;
      coefficients[2] = 0;
      coefficients[3] = 0;
      return;
    }
    if ((p0.y == p1.y) && (p0.z == p1.z) && (p0.x != p1.x)) {
      // y and z have constant values
      // line is defined by y = B, z = C
      line_flag = Flag.YZ_CONSTANT;
      coefficients[0] = p0.y;
      coefficients[1] = p0.z;
      coefficients[2] = 0;
      coefficients[3] = 0;
      return;
    }
    if ((p0.x == p1.x) && (p0.y != p1.y) && (p0.z != p1.z)) {
      // x has a constant value
      // plane (degenerate) is defined by x = A, z = By + C
      line_flag = Flag.X_CONSTANT;
      coefficients[0] = p0.x;
      coefficients[1] = (p1.z-p0.z)/(double)(p1.y-p0.y);
      coefficients[2] = (p0.z*p1.y-p0.y*p1.z)/(double)(p1.y-p0.y);
      coefficients[3] = 0;
      return;
    }
    if ((p0.y == p1.y) && (p0.x != p1.x) && (p0.z != p1.z)) {
      // y has a constant value
      // plane (degenerate) is defined by y = B, z = Ax + C
      line_flag = Flag.Y_CONSTANT;
      coefficients[0] = p0.y;
      coefficients[1] = (p1.z-p0.z)/(double)(p1.x-p0.x);
      coefficients[2] = (p0.z*p1.x-p0.x*p1.z)/(double)(p1.x-p0.x);
      coefficients[3] = 0;
      return;
    }
    if ((p0.z == p1.z) && (p0.x != p1.x) && (p0.y != p1.y)) {
      // z has a constant value
      // plane (degenerate) is defined by z = C, y = Ax + B
      line_flag = Flag.Z_CONSTANT;
      coefficients[0] = p0.z;
      coefficients[1] = (p1.y-p0.y)/(double)(p1.x-p0.x);
      coefficients[2] = (p0.y*p1.x-p0.x*p1.y)/(double)(p1.x-p0.x);
      coefficients[3] = 0;
      return;
    }
    if ((p0.x != p1.x) && (p0.y != p1.y) && (p0.z != p1.z)) {
      // no variables have a constant value
      // plane (degenerate) is defined by y = Ax + B, z = Cx + D
      line_flag = Flag.NO_CONSTANT_VARIABLES;
      coefficients[0] = (p1.y-p0.y)/(double)(p1.x-p0.x);
      coefficients[1] = (p0.y*p1.x-p0.x*p1.y)/(double)(p1.x-p0.x);
      coefficients[2] = (p1.z-p0.z)/(double)(p1.x-p0.x);
      coefficients[3] = (p0.z*p1.x-p0.x*p1.z)/(double)(p1.x-p0.x);
    }
  }

  /**
   *  Calculates new coefficients that for the new point.  Uses the
   *  new coefficients if they are relatively close to to the previous
   *  ones.  Kills off the invariant if they are not.
   *
   *  @return true if the new equation worked, false otherwise.
   **/
  public boolean try_new_equation (long x, long y, long z) {

    // Calculate max separation using this point and the existing 3 points.
    def_points[3] = new Point (x, y, z);
    double sep = maxsep_triples (def_points);

    // If this point increased the separation, recalculate a, b, and c.
    if (sep > separation) {
      separation = sep;
      double[] coef;
      try {
        // Requires that the first 3 elements of def_points must be non-null.
        coef = calc_tri_linear (def_points);
        if (Debug.logDetail())
          wrapper.log ("Calc new plane with points " + def_points[0] + " " +
                       def_points[1] + " " + def_points[2] + " " +
                       def_points[3]);
      } catch (Exception e) {
        return (false);
      }

      // if the a, b, or c is a new min/max remember it.
      if (coef[0] < min_a) min_a = coef[0];
      if (coef[0] > max_a) max_a = coef[0];
      if (coef[1] < min_b) min_b = coef[1];
      if (coef[1] > max_b) max_b = coef[1];
      if (coef[2] < min_c) min_c = coef[2];
      if (coef[2] > max_c) max_c = coef[2];
      if (coef[3] < min_d) min_d = coef[3];
      if (coef[3] > max_d) max_d = coef[3];

      // Pick a new a, b, and c as the average of their endpoints
      a = (min_a + max_a) / 2;
      b = (min_b + max_b) / 2;
      c = (min_c + max_c) / 2;
      d = (min_d + max_d) / 2;
      if (Invariant.logOn() || debug.isLoggable(Level.FINE))
        wrapper.log (debug, wrapper.ppt.name() + ": Trying new a (" + a +
                     "), b (" + b + "), c (" + c + "), and d (" + d + ")");
      wrapper.log ("min max: Trying new a (" + a +
                   "), b (" + b + "), c (" + c + "), and d (" + d + ")");
      // if the new coefficients are 'equal' to their min and max and
      // this point fits, then this new equation is good enough both
      // for existing points and the new point.
      if (Global.fuzzy.eq(a, min_a) && Global.fuzzy.eq(a, max_a) &&
          Global.fuzzy.eq(b, min_b) && Global.fuzzy.eq(b, max_b) &&
          Global.fuzzy.eq(c, min_c) && Global.fuzzy.eq(c, max_c) &&
          Global.fuzzy.eq(d, min_d) && Global.fuzzy.eq(d, max_d) &&
          (Global.fuzzy.eq (-d, a*x+b*y+c*z))) {
        if (debug.isLoggable(Level.FINE))
          debug.fine (wrapper.ppt.name() + ": New a (" + a + ") and b ("
                      + b + ") and c (" + c + ")");
        return (true);
      } else {
        return (false);
      }
    } else { // this point doesn't increase the separation

      return (false);
    }
  }

  /**
   * Calculates the separation between p1 and p2.
   *
   * @param p1  First point
   * @param p2  Second point
   *
   * @return the distance between p1 and p2, or 0 if either point is null,
   * or NaN if either point contains a NaN
   */
  double separation(Point p1, Point p2) {

      // These variable Types are double so the values won't wrap around.

    // make sure both points are specified
    if ((p1 == null) || (p2 == null)) {
      return (0);
    }

    double xsep = (p1.x - p2.x);
    double ysep = (p1.y - p2.y);
    double zsep = (p1.z - p2.z);
    return Math.sqrt(xsep*xsep + ysep*ysep + zsep*zsep);
  }

  /**
   * Calculates the three points that have the maximum separation in pa and
   * places them as the first three elements of pa.
   *
   * @param pa  Array of points.  Must have at least 3 elements.  Can be any
   *            length and can contain nulls (which will be ignored).
   *            Is side-effected so that the first three elements contain the
   *            points with the maximum total separation; this may introduce
   *            duplicates into the array.
   *            The first three elements are null if no triple of elements
   *            is NaN-free.
   *
   * @return the maximum separation found.
   */
 double maxsep_triples (/*@Nullable*/ Point[] pa) {

   // cache values for the (square of the) distance between each pair of
   // points, to avoid duplicating work
   double[][] separations = new double[pa.length][pa.length];
   for (int i=0; i<pa.length-1; i++) {
     for (int j=i+1; j<pa.length; j++) {
       separations[i][j] = separation(pa[i], pa[j]);
     }
   }

   Point p1 = null, p2 = null, p3 = null;

   // max_separation is the separation metric for the most separated
   // triple of points.  We use the sum of the separations as the
   // metric.  (The metric "min of the three separations" does not work
   // because it doesn't choose a unique set of points, making the
   // result dependent on the order in which the points are seen; more
   // seriously, it may choose a collinear set of points even when a
   // non-collinear set exists.)
   double max_separation = Double.MIN_VALUE;

   for (int i=0; i < pa.length-2; i++) {
     for (int j=i+1; j < pa.length-1; j++) {
       double sep_i_j = separations[i][j];
       for (int k = j+1; k < pa.length; k++) {

         // using the sum of separations is valid, as long as
         // separation is the actual distance between points and not
         // the square of the distance
         double separation = sep_i_j + separations[i][k] + separations[j][k];
         if (separation > max_separation) {
           max_separation = separation;
           p1 = pa[i];
           p2 = pa[j];
           p3 = pa[k];
         }
       }
     }
   }

   pa[0] = p1;
   pa[1] = p2;
   pa[2] = p3;

   if (Debug.logDetail())
     wrapper.log ("maxsep_triples = " + pa[0] + " " + pa[1] + " " + pa[2]);
   return (max_separation);
 }

  // Given ((x0,y0,z0),(x1,y1,z1), (x2,y2,z2), calculate a, b, c and d
  // such that ax + by + cz + d = 0.
  // a, b, c and d are mutually prime, integers and a is positive.
  //
  // A visual explanation of the math can be found in:
  // http://www.efm.leeds.ac.uk/CIVE/CIVE2599/cive2599-summary-overheads-full.pdf
  //
  // The standard form of a plane can be calculated using a vector
  // normal "n" <n1, n2, n3> to the plane and a vector from the origin
  // to the plane "r" e.g. <x0, y0, z0>.
  //
  // The vector normal to the plane can be calculated by doing the cross product
  // of two line segments on the plane (<point2 - point0> and <point1 - point0>)
  //
  // Given the normal vector and a point on the plane, the standard form
  // of the plane is:
  // n1*x + n2*y + n3*z = p where p is the dot product of r and n
  //
  //

  /**
   * Calculates the coefficients (a, b, c) and constant (d)
   * for the equation ax + by + cz + d = 0 using the
   * first three points in points.
   *
   * @param points  array of points to use to calculate the coefficents.  Only
   *                the first three points (all of which must be non-null) are used.
   *
   * @return a four element array where a is the first element, b the
   * second, c the third, and d is the fourth.  All elements are
   * mutually prime, integers and a is positive
   */
  // TODO: should just pass in the first three elements rather than passing
  // in the point array.
  public double[] calc_tri_linear(/*@Nullable*/ Point[] points) { // TODO: remove annotation after refactoring method

    Point p0 = points[0];
    Point p1 = points[1];
    Point p2 = points[2];

    // the following code is taken from the source page of
    // http://jblanco_60.tripod.com/plane.html (the math behind the
    // implementation is detailed above)

    long px = p0.x;
    long py = p0.y;
    long pz = p0.z;
    long qx = p1.x;
    long qy = p1.y;
    long qz = p1.z;
    long rx = p2.x;
    long ry = p2.y;
    long rz = p2.z;

    // define the i, j, k components of the two line segments on the plane
    long a1 = (qx - px);
    long a2 = (qy - py);
    long a3 = (qz - pz);

    long b1 = (rx - px);
    long b2 = (ry - py);
    long b3 = (rz - pz);

    // calculate the i, j, k components of the cross product
    long i_var = (a2 * b3) - (b2 * a3);
    long j_var = -((a1 * b3) - (b1 * a3));
    long k_var = (a1 * b2) - (b1 * a2);

    // Calculate the value of the constant.  Note that given the format
    // of the point-normal form, we multiply by the negative of P(x,y,z)

    long Stand_const = ((i_var*(-px)) + (j_var*(-py)) + (k_var*(-pz)));

    // Put the coefficients in lowest terms by dividing by the gcd.
    // If the gcd is 1, no harm done -- divide by 1.
    // Type is double to ensure floating-point division.
    double myGCD = MathMDE.gcd(MathMDE.gcd(i_var, j_var) , MathMDE.gcd(k_var, Stand_const));
    double standard_i_var = i_var/myGCD;
    double standard_j_var = j_var/myGCD;
    double standard_k_var = k_var/myGCD;
    double standard_Stand_const = Stand_const/myGCD;

    if (wrapper != null) {
      wrapper.log("Preliminary: " + i_var + " " + j_var + " GCD: " + MathMDE.gcd(i_var, j_var));
      wrapper.log("Preliminary: " + k_var + " " + Stand_const + " GCD: " + MathMDE.gcd(k_var, Stand_const));
      wrapper.log("GCD: " + myGCD);
    }

    // Ensure postive x coefficient by negating if "x" term is negative.
    if (i_var < 0) {
      standard_i_var = - standard_i_var;
      standard_j_var = - standard_j_var;
      standard_k_var = - standard_k_var;
      standard_Stand_const = - standard_Stand_const;
    }

    double[] coef = new double[] {standard_i_var,
                                  standard_j_var,
                                  standard_k_var, standard_Stand_const };

    return coef;
  }

  public boolean enoughSamples() {
    return isActive();
  }

  // If the values form a line,
  // we have no confidence that it's a plane.
  // The line is less interesting because the invariant
  // is already defined by a linear binary invariant.
  public double computeConfidence() {
    if (isActive()) {
      return Invariant.conf_is_ge(values_seen, MINTRIPLES);
    } else
      return 0;
  }

  public String repr() {
    return "LinearTernaryCore" + wrapper.varNames() + ": "
      + "a=" + a
      + ",b=" + b
      + ",c=" + c
      + ",d=" + d
      + ",values_seen=" + values_seen;
  }

  public String point_repr(Point p) {
    if (p == null)
      return "null";
    else
      return "<" + p.x + "," + p.y + "," + p.z + ">";
  }

  public String cache_repr() {
    StringBuffer result = new StringBuffer();
    for (int i=0; i<MINTRIPLES; i++) {
      if (i!=0) result.append("; ");
      result.append(point_repr(def_points[i]));
    }
    return result.toString();
  }

  // In this class for convenience (avoid prefixing "LinearBinaryCore").
  static String formatTerm(double coeff, /*@Nullable*/ String varname, boolean first) {
    return LinearBinaryCore.formatTerm (coeff, varname, first);
  }

  public String format_using(OutputFormat format,
                             String vix, String viy, String viz,
                             double a, double b, double c, double d) {

    if ((a == 0) && (b == 0) && (c == 0) && (d == 0)) {
      return wrapper.format_too_few_samples(format, null);
    }

    if (format == OutputFormat.SIMPLIFY) {
      return format_simplify(vix, viy, viz, a, b, c, d);
    }

    if ((format.isJavaFamily()))
      {

        return formatTerm(a, vix, true)
          + formatTerm(b, viy, (a == 0))
          + formatTerm(c, viz, (a == 0) && (b == 0))
          + formatTerm(d, null, (a == 0) && (b == 0) && (c == 0))
          + " == 0";

      }

    if ((format == OutputFormat.DAIKON)
        || (format == OutputFormat.ESCJAVA))
      {
        String eq = " == ";
        return formatTerm(a, vix, true)
          + formatTerm(b, viy, (a == 0))
          + formatTerm(c, viz, (a == 0) && (b == 0))
          + formatTerm(d, null, (a == 0) && (b == 0) && (c == 0))
          + eq + "0";
      }

    throw new Error("unrecognized output format " + format);
    // return null;
  }

  public static String format_simplify (String vix, String viy, String viz,
                                       double da, double db, double dc, double dd) {
    int ia = (int) da;
    int ib = (int) db;
    int ic = (int) dc;
    int id = (int) dd;

    String str_a, str_b, str_c, str_d;
    if (ia != da || ib != db || ic != dc || id != dd) {
      // floating point

      // Disabled for the moment, since mixing integers and floating
      // literals seems to give Simplify indigestion. For instance:
      // (BG_PUSH (<= w 3))
      // (BG_PUSH (EQ 0 (* w 2.0d0)))
      // (<= w 1)

      // str_a = Invariant.simplify_format_double(da);
      // str_b = Invariant.simplify_format_double(db);
      // str_c = Invariant.simplify_format_double(dc);
      // str_d = Invariant.simplify_format_double(dd);
      return "(AND)"; // really s/b format_unimplemented
    } else {
      // integer
      str_a = Invariant.simplify_format_long(ia);
      str_b = Invariant.simplify_format_long(ib);
      str_c = Invariant.simplify_format_long(ic);
      str_d = Invariant.simplify_format_long(id);
    }

    String str_z = viz;
    String str_x = vix;
    String str_y = viy;
    String str_ax = (da == 1.0) ? str_x : "(* " + str_a + " " + str_x + ")";
    String str_by = (db == 1.0) ? str_y : "(* " + str_b + " " + str_y + ")";
    String str_cz = (dc == 1.0) ? str_z : "(* " + str_c + " " + str_z + ")";
    String str_axPbyPcz = "(+ " + str_ax + " " + str_by + " " + str_cz + ")";
    String str_axPbyPczPd = (dd == 0.0) ? str_axPbyPcz :
      "(+ " + str_axPbyPcz + " " + str_d + ")";
    return "(EQ 0 " + str_axPbyPczPd + ")";
  }

  public String format_using(OutputFormat format,
                             String vix, String viy, String viz
                             ) {
    String result = format_using(format, vix, viy, viz, a, b, c, d);
    if (result != null) {
      return result;
    }

    return wrapper.format_unimplemented(format);
  }

  // // Format as "x = cy+d" instead of as "y = ax+b".
  // public String format_reversed(String x, String y) {
  //   assert a == 1 || a == -1;
  //   return format(y, x, a, -b/a);
  // }

  public boolean isSameFormula(LinearTernaryCore other) {

    // If both have yet to see enough values
    if (!isActive() && !other.isActive()) {

      // Same formula if all of the points match
      if (values_seen != other.values_seen)
        return (false);
      for (int ii = 0; ii < values_seen; ii++) {
        if (def_points[ii] != other.def_points[ii])
          return (false);
      }
      return (true);
    } else {
      return ((values_seen >= MINTRIPLES)
              && (other.values_seen >= MINTRIPLES)
              && (a == other.a)
              && (b == other.b)
              && (c == other.c)
              && (d == other.d));
    }
  }

  public boolean isExclusiveFormula(LinearTernaryCore other) {
    if (!isActive() ||
        !other.isActive()) {
      return false;
    }

    return ((a == other.a)
            && (b != other.b)
            && (c != other.c)
            && (d != other.d));
  }

  /**
   * In general, we can't merge formulas, but we can merge invariants with
   * too few samples to have formed a plane with invariants with enough
   * samples.  And those will appear to have different formulas.
   */
  public boolean mergeFormulasOk() {
    return (true);
  }

  /**
   * Merge the linear ternary cores in cores to form a new core. Any
   * core in the list that has seen enough points to define a
   * plane, must define the same plane. Any cores that have not
   * yet seen enough points, will have each of their points applied to
   * that invariant.  The merged core is returned.  Null is
   * returned if the cores don't describe the same plane
   *
   * @param cores   List of LinearTernary cores to merge.  They should
   *                all be permuted to match the variable order in
   *                ppt.
   */
  public /*@Nullable*/ LinearTernaryCore merge (List<LinearTernaryCore> cores, Invariant wrapper) {

    // Look for any active planes.  All must define the same plane
    LinearTernaryCore first = null;
    for (LinearTernaryCore c : cores) {
      if (!c.isActive())
        continue;
      if (first == null)
        first = c.clone();
      else {
        if (!Global.fuzzy.eq (first.a, c.a)
            || !Global.fuzzy.eq (first.b, c.b)
            || !Global.fuzzy.eq (first.c, c.c)
            || !Global.fuzzy.eq (first.d, c.d))
          return (null);
      }
    }

    // If no active planes were found, created an empty core
    if (first == null)
      first = new LinearTernaryCore (wrapper);
    else
      first.wrapper = wrapper;

    // Merge in any points from non-active cores
    // Note that while null points can't show up in a core that is not active
    // because it hasn't seen enough values, they can show up if the core
    // is not active because it is a line.
    for (LinearTernaryCore c : cores) {
      if (c.isActive())
        continue;
      for (int j = 0; j < c.values_seen; j++) {
        Point cp = c.def_points[j];
        if (cp == null)
          continue;
        if (Debug.logDetail())
          wrapper.log ("Adding point " + cp + " from " + c.wrapper.ppt);
        InvariantStatus stat = first.add_modified (cp.x, cp.y, cp.z, 1);
        if (stat == InvariantStatus.FALSIFIED)
          return (null);
        //if (wrapper.is_false())
        //  return (null);
      }
    }

    return (first);
  }

}
